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Abstract Wc present a series of numerical simulations of the quiet Sun plasma 
threaded by magnetic fields that extend from the upper convection zone into the 
low corona. We discuss an efficient, simplified approximation to the physics of 
optically thick radiative transport through the surface layers, and investigate 
the effects of convective turbulence on the magnetic structure of the Sun's 
atmosphere in an initially unipolar (open field) region. We find that the net 
Poynting flux below the surface is on average directed toward the interior, while 
in the photosphere and chromosphere the net flow of electromagnetic energy is 
outward into the solar corona. Overturning convective motions between these 
layers driven by rapid radiative cooling appears to be the source of energy for 
the oppositely directed fiuxes of electromagnetic energy. 

Keywords: Convection; Corona; Magnetic fields; Photosphere; Radiative Trans- 
fer 



1. Introduction 

To understand the physics of solar activity, we must understand the magnetic 
and energetic connection between the Sun's convective envelope and corona. 
The magnetic fields that mediate or energize most, if not all, solar activity are 
generated below the visible surface within the turbulent convection zone. Yet 
most of what we can directly measure originates from the solar atmosphere, 
where physical conditions are fundamentally different from that of the interior. 
While helioseismic inversions provide an invaluable window into the physics of 
the Sun's interior, understanding the physical connection between subsurface 
features and those observed in the solar atmosphere requires a realistic forward 
model. 

But what level of realism in a numerical model is necessary to describe the 
complex magnetic connectivity and energetics of the solar atmosphere lying 
between the visible surface and the corona? It is of great benefit, for example, 
to formulate a simple, well-defined problem, and set up an idealized numerical 
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experiment that sheds hght on the relevant physical processes in an otherwise 
complex system. In this way, important progress has been made in our under- 
standing of the physics of magnetic flux emergence in highly stratified model 
atmospheres (see, e.g., [Manchester et al\ 120041 [Murray et al.\ 120061 |Magara[ 
[Galsgaard et a/.[ [MTTl iFaiH [MM lArchontis and Hoodl [^imi|) . 

Yet the observed evolution of the photospheric magnetic field is often far more 
complex, particularly in and around CME- and flare-producing active regions. 
It is difficult to set up a simple magnetic and energetic configuration and an 
associated physics-based photospheric boundary condition that can initialize a 
simulation of the solar atmosphere and faithfully mimic the coronal evolution 
of a complex active region. If we wish to perform first-principles quantitative 
studies of phenomena such as eruptive events, the energization of the solar 
wind, active region decay, the transport of magnetic free energy and helicity 
into the solar atmosphere, and the physics of coronal heating, it is essential to 
evolve a turbulent model convection zone and corona within a single, large-scale 
computational domain. 

To achieve this, we must accommodate the fundamental energetics of the 
system while still retaining the ability to study the interplay between large and 
small-scale magnetic structures that evolve over different timescales. Clearly, 
radiative transport plays a critical role in the energy balance of the atmospheric 
layers that bridge the gap between the visible surface and corona. Surface cool- 
ing drives convection, and convective turbulence both generates magnetic field 
and mediates the flux of magnetic energy that enters the solar atmosphere. 
Yet the physics of radiative transport can be computationally expensive to 
treat realistically, even in the context of small-scale domains that do not in- 
clude the convection zone and corona within a single computational volume. 
For example, energetically important transitions in the solar chromosphere are 
often decoupled from the local thermodynamic state of the plasma (a state 
of non-local thermodynamic equilibrium, or non-LTE) suggesting that a truly 
realistic numerical model must also couple the macroscopic radiative transfer 
and level population equations to the system of conservation equations (see, 
e.g., |McClymont"aH d Canficld, 1983] [Fisher, C anficld , and McClymont[ [IM51 



ICarlsson and Stein[ri992i|Abbett and Hawley[ll999llAlIred~eran[MJ5)) . To com- 
plicate matters further, non-thermal physics, and the physics of ion-neutral drag 
may substantially affect the energy balance of the chromosphere ( [Krasnoselskikh et al., 2010] ) . 

While it remains impractical to perform large-scale, 3D, non-LTE radia- 
tive MHD calculations without employing substantial approximations to make 
the system tractable, it is now common practice to realistically treat optically 
thick surface cooling in the upper convection zone and photosphere in LTE (a 
good approximation in these layers). There are many examples of thin-layer, 
high-resolution calculations that incorporate solutions to the non-gray radiative 
transfer equation in Cartesian domains that include the upper convection zone 
and extend into the low chromosphere ([Bercikl 120021 IStein and Nordlundl 120061 



Georgobiani et a/."] 120071 [Rempel, Schiissler, and Knolker[ 120091 [Cheung et al. 



2010p . In addition, calculations that realistically treat radiative transfer have 



been applied to simulations of solar granulation in relatively small-scale domains 

that also include a transition region and corona ( [Martinez-Sykora, Hansteen, and Carlsson] 
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120081 [Martmez-Sykora, Hansteen, and Car lsson[ 120091 |Carlsson, Hansteen, and Gudiksenl 

mm . 

Our goal, however, is to expand the size of such computational domains to 
active region or even global spatial scales while still retaining as realistic a ther- 
modynamic environment as is feasible. Thus, we strive to develop the simplest 
model possible that allows us to capture the essential physics of the convection- 
zone-to-corona system while still maintaining the computational efficiency of 
models in which optically thick radiative cooling is treated in a parameterized 
fashion (e.g., lAbbetti [20071 [Fang et al.\ I2010ap . In this way, we hope to make 



practical the performance of physics-based, first principles simulations, allowing 
for quantitative, parameter-space studies of processes such as filament formation, 
active region emergence and decay, and flare and CME initiation. 

To simultaneously evolve a realistic model convection zone and corona at any 
spatial scale presents a number of daunting challenges. The upper convection 
zone and low solar atmosphere are highly stratified — average thermodynamic 
quantities change by many orders of magnitude as the domain transitions from a 
relatively cool, turbulent regime below the visible surface, to a hot, magnetically- 
dominated and shock-dominated regime high in the corona. The physics of the 
gas transitions from a high-/? plasma where the magnetic field is advected by 
the gas (away from strong active region complexes) to a low-/? regime where the 
gas is constrained to move along magnetic field lines. In addition, the radiation 
field transitions from being optically thick to optically thin. Temporal and spatial 
scales are highly disparate. Large concentrations of magnetic flux are compressed 
within intergranular lanes and evolve at convective turnover timescales, while 
large coronal loops form and persist for days as active regions emerge and evolve 
over a course of many months. In addition, the large-scale magnetic structure of 
the corona can change in a fraction of a second, as small-scale localized magnetic 
reconnection suddenly reorganizes the large-scale fleld, often triggering eruptive 
events along the way. 

The corona presents particular challenges. It is well known that in order to 
accurately reflect the thermodynamics of this region, a model should include 
the effects of electron heat conduction along magnetic field lines and radiative 
cooling in the optically thin "coronal approximation" . In addition, some physics- 
based {e.g., Joule heating) or empirically based source of coronal heating must 
be present (often introduced at the lower photospheric boundary) if the model 
corona is to remain hot. But to generate a realistic magnetic carpet, and to study 
the interaction of granular convection with coronal structures, requires there to 
be a turbulent model convection zone, and therefore some form of optically thick 
surface cooling. 

In lAbbettl (|2007p we introduced this physics into a 3D MHD convection-zone- 
to-corona model in the simplest, most computationally efficient way possible — 
we simply ignored the optically thick radiative transfer equation entirely, and in- 
stead used a parameterized Newton cooling function carefully calibrated against 
smaller-scale, more realistic radiative-MHD models of magneto-convection where 
the frequency-dependent LTE transfer equation was solved along with the MHD 
system ( |Bercik, 2002] ) . 
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This approach has been successful in studying the structure of quiet-Sun 
magnetic fields and active region flux emergence (jAbbetti 120071 |Fang et all 
I2010al [Fang et aT]l2010b| . Yet this treatment, while computationally efficient, 
has a number of limitations. Its principle drawback is that it is ultimately ad 
hoc and requires other, more realistic simulations as a basis for calibration in 
order to get meaningful results. The simplified cooling is imposed at a particular 
height or over a range of gas density, and is not generated in a physical way 
as a function of optical depth. To address these limitations, we build upon a 
technique introduced by lAbbett and Fisherl (j2010p , and in Section [3] derive a 
simple, flux-conservative approximation to optically thick cooling that is based 
on the gray radiative transfer equation in LTE. We then incorporate a form of 
this efficient, physics-based approximation into the RADMHD convection-zone- 
to-corona model of lAbbettl l|2007p . which we briefly describe in Section [21 In 
Section [2 we present new models of an open-field coronal hole region, and study 
the transport of magnetic energy from below the surface into the corona. Finally, 
in Section [5] we summarize our results. 



2. Numerical Methodology 

The parallel code RADMHD solves the following MHD conservation equations 
semi-implicitly on a three-dimensional Cartesian mesh: 



dp 

dt 



dpu 



V- (pu) = 0, (1) 

Pg, (2) 



/ BB 

— + V-(uB-Bu) = -Vx(ryVxB), (3) 

-^ + V-(eu) = -pV-u+-^|VxB|2 + $ + Q,. (4) 
ot 47r 

The components of the state vector have the usual definitions: p, u, e, p, 
B, and g denote the gas density, velocity, internal energy per unit volume, gas 
pressure, magnetic field, and gravitational acceleration respectively. Here, we 
assume Gaussian units. The viscous stress tensor is assumed to be of the form 
Ilij = 2pv[Dij — 1/3(V • u)5ij], where Dij = \/2{dui/dxj + duj/dxi) and 5ij 
denotes the Kronecker delta function. The function 4> = ■ IlijDij represents 
the rate of energy dissipation through viscous diffusion, and ly and r/ refer to the 
coefficients of kinematic viscosity and magnetic diffusivity, respectively. These 
coefficients are assumed constant, and are set to values that correspond to the 
grid-scale viscous and resistive dissipation. The source term Q includes impor- 
tant energy sources and sinks such as radiative cooling, the divergence of the 
electron heat flux (in the portion of the domain representing the model transition 
region and corona) , and any desired empirically based coronal heating function. 
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A complete discussion of the components of this energy source term is provided 
bv lAbbettKlMTT)) . The system is closed with a non- ideal equation of state, using 
tabular data provided by the OPAL project ( [Rogers, 2000[ ). In this article, the 
portion of the domain corresponding to the corona is heated by the empirically 
based coronal heating function described in lAbbettl (|2007p . and the effects of 
Joule dissipation within this region are ignored. 

The semi-implicit numerical scheme is parallelized on a domain-decomposed 
mesh, and the core technique is based on operator splitting with a high-order 
Crank-Nicholson temporal discretization. We treat the electron thermal con- 
duction, viscous and Joule dissipation, and radiative losses implicitly using a 
Jacobian-Free Newton-Krylov (JFNK) solver, and require that the remainder 
of the system be treated explicitly using the Central Weighted Essentially Non- 
Oscillatory (CWENO) method of ( [Kurganov and Levy[ [20001 IBalbas and Tadmofl 
I2006p . In this way, we remain Courant limited by the magnetosonic wavespeed, 
and can follow the dynamics of the system in a reliable way (we may choose 
to relax this constraint when evolving active region magnetic fields over longer 
timescales). Any local divergence error introduced into the magnetic field as 
a result of the CWENO central scheme is dissipated by adding an additional 
artificial source term proportional to V(V • B) to the induction equation. A 
detailed description of the numerical methodology employed by RADMHD can 
be found in Section 2 of lAbbettl ([W7| . 

A number of enhancements and improvements have been incorporated into 
the RADMHD source code since its initial release in 2007. Most improvements 
arc in the form of improved performance and robustness, better MPI load bal- 
ancing and scaling, and other enhancements in the code's speed and efficiency. 
Among the enhancements arc: i) a simplified and improved table inversion and 
interpolation algorithm that is necessary to incorporate the OPAL data into the 
code's non-ideal equation of state and the CHIANTI data ( [Young et ai, 2003| ) 
into the code's treatment of optically thin radiative cooling; ii) a new adaptive 
error algorithm in the GMRES (Generalized Minimum RESidual) substep of the 
JFNK solver that greatly improves convergence rates; Hi) a more robust, global 
non-linear CWENO weighting scheme in the explicit substep of RADMHD; and 
iv) an option to evolve log p rather that p itself via the following rewrite of 
Equation ([T|): 

r\ 1 

^-^ + V-(ulnp) = (Inp- l)V-u. (5) 
ot 

Since the model atmosphere is highly stratified, this is often useful as a means 
of making the code more robust, while at the same time retaining the desired 
shock-capture characteristics of the numerical scheme. More details on these and 
other algorithmic improvements will be provided in a technical document under 
preparation for inclusion with the next release of the code. 

In essence, however, the core numerical methods of RADMHD remain the 
same as that presented in lAbbetTl (|2007p . In this article, we focus on the portion 
of the energy source term Q of Equation ^ that contains the approximation 
for optically thick radiative cooling. 
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3. An Approximate Treatment of Optically Thick Cooling 

Radiative cooling drives surface convection and is a crucial contributor to the 
energy balance in the region of the solar atmosphere bridging the convection 
zone and corona. Yet a full frequency-dependent solution to the LTE radiative 
transfer equation can be computationally expensive for large-scale convcction- 
zone-to-corona calculations, particularly for active region or filament models 
where timescales are such that the radiative cooling must be updated at intervals 
close to the MHD CFL limit. Here, we build upon the approach introduced 
bv lAbbett and Fisherl (|2010p . and derive an approximate, frequency- integrated 
expression for optically thick radiative cooling that is based on the gray transfer 
equation in LTE. We begin by considering the net cooling rate for a volume of 
plasma at a particular location in the solar atmosphere: 

R^JdnJd,^{7],~K,I,). (6) 

Here, represents solid angle, and v the frequency. The subscript v indicates 
that the emissivity, opacity, and specific intensity (rji,, k^, and I^, respectively) 
depend on frequency If wc define the source function Si, as the ratio of the 
emissivity to opacity, and rearrange the order of integration, we can recast the 
net cooling rate in the following form 

R= jdvK„jAn{S,-h). (7) 

We define the mean intensity as = (l/47r) jdVH^, and note that the source 
function is independent of direction. This allows Equation ([7]) to be expressed 
as 

R = 4tt J diy (S, - J,) . (8) 

If we now assume a locally plane-parallel geometry, the formal solution for the 
specific intensity can be written as (e.g.. IMihalasi IT978| 

/,(T„Ai)=/ dr' 5,(t'), (9) 

Jo IMI 

where fj. refers to the cosine angle and to the frequency-dependent optical 
depth. Wc can now recast the expression for the mean intensity in terms of an 
integral over optical depth and cosine angle, 

1 r°° /"i p-h-^-T'\/\tJ.\ 

Mru)^T^ dr'SAr) dU,\ , (10) 

^ Jo Jo IMI 

This allows the integral over fi to be evaluated and expressed in terms of an 
exponential integral function, 

1 r°° 

■Mr.) = / dr'SAT')Eii\T, - r'l). (11) 



2 7o 
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Up to now, no approximation other than an assumption of a locaUy plane par- 
allel geometry has been made. We now follow the analysis of lAbbett and Fisherl 
(|2010p and note that the first exponential integral function [Ei{\ti, — r'|)] in 
Equation pip is singular when t' = t^, and that this singularity is intcgrable. 
Since Ei is peaked around t^, contributions from Si,{t') will be centered around 
Sy{T,j). Thus, to lowest order we can approximate the mean intensity by 

1 

JAru) « 7;SAr.) / dT'Ei{\T, - r'|). (12) 

The integral over optical depth is now easily evaluated, and the result is ex- 
pressed in terms of the second exponential integral function [£'2]: 

JM^sM^i'^^y (13) 

Wc now rearrange the terms in the above equation, and substitute 1 — {t^ )I {tv) 
« i?2(T,y)/2 into Equation ([H]), to arrive at an approximation for the net cooling 
rate, 

R^2^ jAuK,S,,E2{T,). (14) 

If wc further assume LTE, the source function can be expressed as the Planck 
function [B,j{T)] coupling the cooling rate to the local temperature of the plasma 
[T]: 

RK2n jdvK,B,{T)E2{T,). (15) 

We now integrate Equation psp over frequency. Since E2(t,^) is bounded 
below by zero and above by unity, the integral in Equation ()15p obeys this set 
of inequalities: 



2KaT'^ >'2tt J Auk^BAT)E2{tA > , 



(16) 



where R is the Planck-weighted mean opacity. 

Because of the range of the E2{x) function, we can use this inequality to write 
the integral in Equation (jTSj) in the form R{f) = 2RC{f)aT'^E2{a{f)f), where in 
general, a is a positive, unknown function of mean optical depth f (dr = —RAz), 
and C(f) is an unknown, f-dependent normalization constant. However, since 
we expect a close relationship between the mean optical depth f and the local 
mean opacity k, we therefore make the ansatz that a is a constant, but with an 
unknown value. The expression for R can then be written 

RKi2CKaT-^E2{aT), (17) 

where C now represents a f-independent normalization constant of integration. 

To determine the normalization constant C, we integrate our cooling function 
from zero to infinity in optical depth over an isothermal slab to obtain the 
total radiative flux. The resulting expression must be equal to the known result 
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^tot = crT"', thus requiring C = a. To calibrate a, we compare the cooling rate 
as a function of depth in test models using this approximation against more 
realistic models of magnetoconvection where the frequency-dependent transfer 



equation is solved in detail (Bercik, 2002). We conclude that the best- fit value is 



a = 1 (see Figure 1 of lAbbett and Fisheri I2010p . This implies that optical depth 
in highly stratified atmospheres is dominated by the local opacity. 

With the parameter a specified, we arrive at an approximation for optically 
thick surface cooling, 

R^2kc7T*E2{t). (18) 

This expression can be efficiently evaluated at each iteration of an MHD calcu- 
lation, and we have implemented this volumetric cooling rate as a part of the 
cell-centered source term Q in Equation (|3|) (the cooling rate R being a negative 
heating rate Q). 

It is possible for the computational grid to be of sufficient resolution to resolve 
the local pressure scale heights of a highly stratified model atmosphere while at 
the same time being poorly resolved in optical depth. This has the potential to 
lead to numerical error in the calculation of the local cooling rate such that 
the total radiative flux may not be conserved. We therefore consider a flux 
conservative formulation similar to the constrained transport schemes common 
to many MHD codes (see IStone and Normanl I1992p . 

We begin by defining a frequency-independent discretized, optical depth for 
each iteration where the MHD state variables are updated. Since our simple 
approximation is based on the assumption of a locally plane-parallel geometry, 
and we are neglecting (for now) the effects of sideways transport, all that is 
required is an integration along the vertical direction {i.e., in the direction of 
the gravitational acceleration). Our discretized expression takes the form 



Ti,j,k-l/2 



l^i.J,k {zn+l/2 - Zn-1/2) (19) 



Here, the grid coordinates i,j,k are defined at cell centers (consistent with the 
centralized numerical scheme implemented in RADMHD), and the optical depth 
is defined at face centers of the mesh cell's control volume perpendicular to 
the z-direction (we now drop the overbar notation, since the above definition 
makes it clear that r is a frequency-averaged quantity) . We use tabular Planck- 
weighted mean opacities (Kij^k) provided by the opacity project ( |Seaton, 2005[ ). 
The coordinate ktop refers to the first ghost cell of the upper coronal boundary 
of the simulation domain, though in practice it is set to an interior cell bounding 
the portion of the domain that represents the optically thin corona. Either way, 
it is presumed that ^.^^ _i/2 = 0. Note that optical depth increases inward 
into the atmosphere in the opposite sense of the height z, which is defined to 
increase outward from the interior toward the visible surface {i.e., dr = —Kdz). 

The radiative cooling of Equation (jl8p can be expressed in terms of a diver- 
gence of a radiative flux. Our treatment of radiative transfer assumes a locally 
plane-parallel geometry, thus we need only consider the radiative flux at the 
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faces of control volumes normal to the vertical direction. This implies that any 
horizontal divergence of the radiative flux is assumed negligible when compared 
to gradients in the vertical direction. The physical justification for this simpli- 
fication is that changes in emissivity and opacity are generally much greater in 
the vertical direction of a highly stratified atmosphere than those expected in 
the transverse direction. This assumption will likely not be valid at the edges 
of sunspots where the lateral emissivity and opacity gradients are expected to 
be large. Given this simplification, the divergence of the radiative fiux can be 
expressed as 

OF 

dz 

We now cast this expression in terms of optical depth dr = —Kdz: 



2R(tT'^E2{t). (20) 



^ = -2aT^E2{r), (21) 

and note that this equation is of the form Q{t) — — /(r)_E2(T) with /(r) = 2<tT'^. 
We now approximate /(r) with a Taylor-series expansion centered about 
accurate up to second order, and reorder the terms so that the expression is of 
the form /(r) = A + Bt: 

fir) - [/(Tfc) - nf'in)] + rfin). (22) 

Here, /(r^) = 2aT^{Tk) refers to the function /(r) evaluated at cell center 
coordinate (i, j, fc), f'{Tk) refers to the function's vertical derivative with respect 
to optical depth evaluated at the same location, and the constants A and B have 
the form A = [/(r/c) — Tfe/'(rfc)] and B = f'{Tk)- For brevity, we have dropped 
the i and j subscripts, but note that these expressions are valid for all grid cells 
at a particular height. 

To obtain the discretized form of Equation ([2T|) , we integrate over the control 
volume of the computational cell. 



F{Tk+i,2)~F{Tk-i,2) = - / '^"\A + BT)E2{T)dr. (23) 

Tfc-l/2 



This integral can be evaluated using the relation di?„(r)/dr = — £'„_i(t) to 
obtain 

FiTk+1/2) - F{Tk-i/2) = A [E:i{Tk+i/ 2) - E:i{Tk- 1/2)] 

+ B [Tk+l/2 E3,{Tk+i/2) - Tk-l/2 E3{Tk_i/2) 

+ Ei{Tk+l/2) - -B4(Tfc_i/2)] . (24) 

All that remains is to define a stencil to evaluate /(t^) and /'(t^) in the expres- 
sions for A and B. RADMHD employs a central scheme, and it is desirable to 
derive a stencil consistent with the formalism of the code. Since the temperature 
T is obtained via a table lookup based on cell-centered values of gas density and 
internal energy per unit volume, we obtain our interpolation stencil by expanding 
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elements of the state vector q{T) in a second-order accurate Taylor series about 

Tfc, 

q{T)=a + b{T-Tk) + ^{T-Tkf, (25) 

and enforce the definition of ccll-avcraged quantities along the z-coordinate axis 
[again, the (i, j) dependence is implicitly assumed, and At = — '''fc-1/2 is 

less than zero], 

1 /■'rfc+1/2 

— / g(r)dr. (26) 

We then expand about each of the points r^-fi, r^, and r^-i; substitute the 
appropriate form of Equation (1251) into Equation (j26p : then perform the integra- 
tion over each respective control volume. This yields a system of equations whose 
solution specifics a, b, and c in terms of known cell averages. The compact stencils 
are equivalent to those of lAbbettl p007|) . and have the form a ~ — (qk+i — 
2qu + gfc-i)/24, b = {qu+i - gfc_i)/(2AT), and c/2 = [qu+i - 2qk + qk-i)/{^T)\ 
The second term of a arises from the integration of the second-order term in the 
Taylor expansion of ^(t), and ensures that the interpolation scheme maintains 
second-order accuracy, and that the following discretized, cell-centered forms of 
A and B have desirable stability properties: 

^ = fk - 7^ ifk+l - 2fk + fk-l) - TkTTT- [fk+l - ./fc-l) (27) 



B=^{fk+i- fk-i). (28) 

Here fk = 2(7 T^J, and = rfe_i/2 — k/j(Az/2). With A and B specified, we 
arrive at an expression for a fiux-conservative approximation to the optically 
thick radiative source term, 

Qinj.k) = A'i^i.,3,k[E3{T, ,j k_i/2) - E3{T, ,j k_^,lf2)] 

+ BKi,j,k [TiJ,k-l/2 Ez[Ti,j.k-l/2) - 'Ti.],k+l/2 £^3 ('''i j,fe+l/2 ) 

+ £^4(Tij,fc_i/2) - £^4(Tij,fc+i/2)] . (29) 

By design, this expression will conserve fiux to machine roundoff, as can easily 
be verified by showing that Qk + Qk~i — ^/c+1/2 ~ -Pfe-3/2 for all points k). 
Thus, we have two ways of implementing our approximation — the flux conser- 
vative approach of Equation (j29|) , and the non-conservative approach obtained 
by directly evaluating Equation (|18p using cell-centered quantities. The flux- 
conserving method requires additional table lookups each iteration to evaluate 
the exponential integrals, but is helpful in cases where the optical depth scale is 
not particularly well-resolved. 

Our approximation for gray LTE cooling is applied only in those regions of the 
computational domain where such an approximation is needed. Specifically, we 
apply the approximation over a range of optical depths that extend from t ~ 10 
to T ~ 0.1. At greater optical depths, wc use the diffusion approximation (with 
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Figure 1. Left: Temperature at the RADMHD model photosphere. Right: Magnetic fieldlines 
threading the low atmosphere over a small sub-domain (the box in the left frame indicates the 
approximate size of the corresponding sub-domain) . The gray slice indicates the average height 
of the visible surface. The domain spans 24 X 24 X 12 Mm^ at a resolution of 512 X 512 X 256. 



tabular Rosseland mean opacities provided by ISeatoni I2005p . and at smaller 
optical depths we use the optically thin approximation (using CHIANTI data 
from IDere et aL[ 119971 and Young et al. 120031 to specify the optically thin cool- 
ing curve) as described in lAbbettl(|2007p and Lundquist, Fisher, and McTiernan 

mm . 

The simulations we present in Section 2] use the non-conservative technique. 
The conservative approach described above was motivated by the fact that in 
some cases, where the model atmosphere is poorly resolved in optical depth, the 
strong cooling prescribed by Equation p8|) can be concentrated in a narrow 
one- or two- zone layer of a model atmosphere. As a practical matter, this 
required that we enforce a limit on the maximum amount of cooling per unit mass 
allowable in any given grid cell. This cooling floor is somewhat artificial, and is 
not necessary in the flux conservative method, which has the effect of spreading 
the cooling over adjoining cells in a more physical way. Another option would 
be to have a separate grid for optical depth, but we decided against this because 
of the possibility of introducing additional interpolation error that is difficult 
to characterize. We are currently testing our new flux-conservative scheme, and 
plan to fully implement it in a new radiation subroutine for RADMHD that also 
includes the important effects of sideways transport. We hope to report on these 
efforts in the near future. 



4. A Model of an Open Flux Region 

We initiate our calculations using the procedure of lAbbettl p007|) . Briefly, we 
begin by relaxing a ID-symmctric average stratification, then expand the domain 
to three dimensions and break the ID symmetry by introducing a small, random 
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energy perturbation in the superadiabatically stratified portion of tlic compu- 
tational domain representing tlie solar convection zone. Convective turbulence 
develops as the simulations progress, and we allow the model convection zone 
to dynamically relax. Wc show results from two separate simulations: one was 
performed locally using 112 processors of a relatively small Beowulf cluster, 
and the other was performed on NASA's Discover supercomputer using 512 
processors. Both simulations simultaneously evolve a model convection zone and 
corona, and each domain has a vertical extent of of 12 Mm, with a 2.5 Mm 
deep model convection zone. The development model that was run on the local 
Beowulf cluster has a domain that spans 21 x 12 x 12 Mm"^ at a resolution of 
448 X 256 X 256, while the larger run on Discover spans 24 x 24 x 12 Mm"^ at a 
relatively high resolution of 512 x 512 x 256. In each case, only « 0.66 percent 
of the total computational effort was expended by the approximate treatment 
of the radiative transfer on average within any given MPI subdomain during a 
given timestep. On the Intel Xeon E5420 CPUs of our local Beowulf cluster, 
the computing time per update of this substep is approximately 0.03 core- 
microseconds per point. The simulations presented here should be considered 
relatively small-scale in the context of the capability of the algorithms presented 
— the code scales well on multiple processors, and once our development work 
is complete we intend to dramatically extend the spatial scale of the models. 

The simulations presented here differ from those of I Abbett I (^0071 in a funda- 
mental way. The approximation we now use for optically thick cooling eliminates 
all of the ad-hoc calibrated parameters present in the older models. Specifically, 
the height and magnitude of the optically thick radiative source term is now 
calculated in a physically self-consistent way based on an optical-depth scale 
rather than on an specified density or height range attenuated by envelope 
functions {c.f. Section 2.1.1 of lAbbetti 12007^ . Once a of Equation has been 
calibrated against more realistic models, no further adjustment is required, and 
each atmosphere relaxes to a state determined by the solution of the system of 
Equations ©-(HI) subject to imposed boundary conditions. 

We apply periodic boundaries in the horizontal directions, and a simple, 
somewhat-artificial closed lower boundary. Specifically, the internal energy per 
unit volume within ghost cells adjacent to the domain's lower boundary is set 
such that a temperature gradient is maintained that best matches the average 
stratification at a corresponding height in the IBercikI (|2002p magnetoconvection 
models. In addition, the ghost cells at the lower boundary are specified such 
that the vertical components of the velocity and magnetic field and the vertical 
gradients of the horizontal components of the velocity and magnetic field are 
zero, and such that the gradient of the gas density is maintained. The upper 
coronal boundary is initially taken to be anti-symmetric during the relaxation 
procedure (i.e., ghost zones are set such that the vertical component of the 
velocity and magnetic field vanishes at the boundary while all other components 
of the MHD state vector maintain a zero vertical gradient across the boundary 
interface) , then is set to a standard zero-gradient boundary condition once mag- 
netic fields are introduced {i.e., ghost cells are set such that all components of the 
MHD state vector maintain a zero vertical gradient across the upper boundary 
interface). For the simulations presented here, once the purely hydrodynamic 
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Figure 2. Temperature in the RADMHD low chromosphere showing a reverse granulation 
pattern. Lighter (darker) colors indicate hotter (cooler) temperatures. In the models, this 
occurs because the radiative cooling diminishes with height, and the pV • v work of converging 
and diverging flows above the intergranular lanes begins to dominate. The horizontal slice spans 
21 X 12 Mm 2 at a resolution of 448 X 256 



model convection zone is relaxed, we introduce a weak 1 G vertically directed 
magnetic field. This is intended to create an open-flux region, such as one might 
expect within a coronal hole. 

As the simulations progress, the convective turbulence acts to stretch and 
amplify the field, and the portion of the domain representing the corona begins 
to heat as a result of the magnetic-ficld-dependcnt empirically-based coronal- 
heating source term. This heating function is based on the IPevtsov et al.\ (I2003P 
power-law relationship between X-ray luminosity and total unsigned magnetic 
flux observed at the surface (see Equation 12 of lAbbettll2007p . For this study, we 
are content to rely on empirical heating rather than Joule dissipation to energize 
the model corona since our focus is on the transport of magnetic energy into the 
atmosphere, not the heating of the low atmosphere and corona. 

Up to the point that magnetic field was introduced into the simulation do- 
main, the model corona was simply an unphysical, cold, nearly evacuated region. 
Once the corona heats, we activate the implicit electron thermal conduction 
source term. This builds a corona, but reduces our timestep somewhat since the 
stiffness of the system is increased and convergence rates in the JFNK substep 
can become an issue. Thus, it tends to be the last step of our relaxation process. 
After several additional turnover times, we begin our analyses. 

The left frame of Figure [T] shows the temperature at the visible surface (t = 
1) of a relaxed convection zone-to-corona model using the new formalism of 
Section [3] The granular pattern and convective turnover times compare well to 
the realistic magnetoconvection models of IBercikI (|2002p . This, along with the 
reverse granulation pattern shown in Figure [21 indicates that our approximation 
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Figure 3. Magnetic fieldlines threading the low atmosphere over a portion of the compu- 
tational domain. The gray slice represents the approximate position of the model's visible 
surface. Horizontally directed magnetic fields due to the spreading of canopy-like structures 
and overturning convective cells permeate the low atmosphere. 



to optically thick radiative transfer is capturing the physics of surface cooling at 
least well enough to generate and sustain solar-like convective features. 

The right frame of Figure [T] shows the complex magnetic structure threading 
a small portion of the simulation domain (as indicated by the cyan box in the 
upper-right corner of the left frame). Figure[3]shows a larger subdomain at a later 
time, and more clearly illustrates the characteristics of the magnetic structure. 
In the region where convective cells turn over, the average plasma-/? remains 
relatively high. As a result, much of the magnetic field remains entrained in the 
plasma, turns over, and is recirculated back below the surface. Thus, at any given 
time, this region is filled with horizontally directed field, and that field tends to 
be less concentrated than is typical of fields entrained in the vortical downdrafts 
present at the visible surface and below. In addition, the presence of canopy-like 
structures (where strong concentrations of field above intcrgranular lanes and 
photospheric downdrafts open into the upper atmosphere and spread out like a 
fan) also contribute to the net amount of horizontally directed magnetic field 
threading the atmosphere below the corona. 

The presence of horizontal fields in the atmosphere has consequences for the 
transport of magnetic energy into the upper atmosphere. Consider the electro- 
magnetic Poynting flux, 

S = — cE X B. (30) 
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Figure 4. The vertical component of the Poynting flux along a layer positioned just below 
optical depth unity (top) and 400 km higher near the tops of overturning granules (bottom). 
Light colors correspond to outward directed flux (toward the corona); dark colors represent 
inward-directed flux (toward the convective interior). Each slice spans 21 X 12 Mm^. 



The vertical component of the Poynting flux is a measure of the amount of 
electromagnetic energy flowing into, or out of the solar atmosphere from below 
the surface where it is generated. In Figure SI we display Sz as a grayscale image 
at two layers in the model atmosphere — dark shades correspond to a flux of 
magnetic energy directed toward the interior, while lighter shades correspond 
to an outward-directed flux. The top frame shows the vertical component of 
Poynting flux along an x-y slice positioned just below the visible surface, and 
the lower frame shows Sz along a slice positioned in the low atmosphere, 400 km 
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Figure 5. The normalized net vertical component of the Poynting flux over a portion of 
the domain centered at the model's photosphere. The dashed vertical line represents the 
approximate height of the visible surface. Above the visible surface, electromagnetic energy 
tends to flow outward toward the corona, while below the surface, energy flows inward toward 
the convective interior. Above z 0.5 (in the model's low chromosphere) the Poynting flux 
tends to remain outwardly directed, but its magnitude is, on average, less than a percent of 
its maximum value. 

higher. Careful examination of Figure S] reveals an imbalance in the outward- 
and inward-directed flux. Below the surface, there appears to be a net flow of 
magnetic energy into the convective interior along the strong vortical downdrafts 
contained within intergranular lanes. Conversely, in the low atmosphere, the ver- 
tical component of the Poynting flux appears more diffuse, and there appears to 
be a net excess of outward directed flux, particularly within overturning granules. 

This is shown more clearly in Figure [5] where Sz is integrated over each layer in 
the computational domain, and plotted as a normalized quantity as a function of 
height \z\ . The dashed vertical line in the figure represents the average height of 
the visible surface. What is clear, is that magnetic energy on average is directed 
downward into the interior below the visible surface. It is at the surface and 
above that the net vertical Poynting flux changes sign and becomes outwardly 
directed. This suggests that the kinetic motion of overturning granules in the 
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model's overshoot layer provides the source of magnetic energy for the corona, 
not the deeper layers below the optical surface, where magnetic flux and energy 
are being pumped down into the interior along intergranular downflows. 

This can be understood in a fairly straightforward way. The vertical compo- 
nent of the Poynting flux can be expressed as 

5, = -^(cE,,xB,,)-z, (31) 
An 

where E/i and refer to the horizontal components of the electric and magnetic 
field respectively. If we assume ideal MHD, then the horizontal component of 
the electric field can be written as follows: 

cEh = —Uyix'Bh — UhxBzZ. (32) 

Just above the visible surface, the magnetic field remains entrained in the fluid 
as convective cells overturn. At this height, the strongest field concentrations are 
located near the edges of overturning granules as divergent flows from neighbor- 
ing cells compress the field. On average, there is more of a contribution to the 
horizontal electric field from the second term of Equation ([5^ , cE^ = — x S^z, 
since the magnetic field becomes more vertical as converging flows compress flux 
into a relatively small area. The contribution of this term to the vertical com- 
ponent of the Poynting flux can be expressed as 4:TrSz = [(— u/iXi3^z)xB/i]- z. 
Simplified, this becomes 4:nSz = —B^CBh- Uh). 

To illustrate the correlation between the horizontal magnetic fields [Bji] and 
the converging surface flows [u/j], consider a weak, vertically oriented, untwisted 
magnetic fiux tube that passes through the surface. Suppose it is acted on by 
a strong converging fiow in a thin layer at the surface. If the magnetic field 
of the tube is oriented in the positive z direction, then just above the surface, 
the compression will tilt the ficldlines and create horizontal components of the 
magnetic field in the opposite direction of the converging fiow. If the magnetic 
field is oriented in the negative z direction, then the horizontal components 
of the field will be aligned with the fiow. Either way, AttSz = —BzCBh - Uh) is 
positive above the surface. Obviously, the dynamics of the model arc far more 
complex than this simple thought experiment. Nevertheless, we do find a net 
positive contribution to the Poynting fiux from the second term of Equation 
([5^ along the edges of overturning granules above the surface where the field is 
being compressed. 

Below the photosphere, the situation is quite different. The strongest magnetic 
fields arc concentrated within localized vortical downdrafts. The asymmetry 
between these strong downdrafts and the broad upwelling plasma in stratified 
convection is well known, and may provide a mechanism whereby magnetic flux 
can be pumped into the interior ([Tobias et al.ll^UUT^ . We find more of a contri- 
bution to the net Poynting flux from the first term of Equation ([32]) . AttSz = 
[(— U2ZxBft)xB/i]- z. This can be expressed more simply as 4:nSz — UzB^, and in 
this form, it is easy to see that a net downward transport of magnetic fiux is con- 
sistent with a net downward-directed Poynting flux. This downward-directed flux 
of electromagnetic energy below the surface is consistent with other simulations 
of radiative-magnetoconvection (see, e.g., |V6gler and Schiissler[ I2007p . 
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The first to recognize this change in direction of tlic fiow of electromagnetic en- 
ergy was lSteiner et a?.1(|2008p who referred to the visible surface as "a separatrix 
for the vertically directed Poynting Flux" . Our results are consistent with their 
findings, although we conclude that in the larger domain, the upward-directed 
net flow of magnetic energy tends to arise from the action of the compressive 
flows of overturning convection as magnetic flux is expelled from cell centers 
and concentrated into the intergranular regions. However, higher in the model 
atmosphere as the gas transitions to a \ow-/3 regime (the upper chromosphere- 
transition region boundary), we also see a small buildup of magnetic flux, for 
reasons similar to those of lSteiner et aLI(|2008p . Namely, that the dynamic chro- 
mosphere transitions to a stable, subadiabatic, magnetically-dominated regime, 
and there is a magnetic reservoir on average as magnetic flux that is advcctcd 
upward enters the stable regime and does not get recirculated back into the 
convective interior (similar in some ways to the overshoot layer at the base of 
the convection zone). This is reflected in the small peak in Figure [5] at a height 
of w 1.3 Mm above the visible surface. Above this transition region interface 
in the open flcld of the model corona, energy is transported via magnetosonic 
and Alfven waves. We note that in these simulations, the magnetic fleld has yet 
to fully saturate (i.e., there is a small increase in magnetic energy over time as 
magnetic field is stretched and amplified by convective turbulence). While this 
indicates that the atmosphere has yet to fully relax, this increase in magnetic 
energy (and any Joule heating below the corona) is negligible in comparison to 
the divergence of the Poynting flux and the work done on the magnetic field by 
convective motions. 

In some sense, the height at which the transition between outward and inward 
flow of electromagnetic energy takes place is less important than the fact that 
such a transition exists. What the simulations seem to suggest is that in quiescent 
regions away from particularly strong concentrations of magnetic flux, there is 
not a continuous flow of electromagnetic energy from below the surface out 
into the corona. Instead, the mechanical energy of convection mediates the flow 
of magnetic energy in the relatively high-/? surface layers where the magnetic 
field remains frozen into the plasma. Of course, in and around very strong 
concentrations of magnetic flux, the situation is undoubtedly quite different. 

5. Discussion and Conclusions 

We have developed an approximate treatment of optically thick radiative surface 
cooling that successfully reproduces the average thermodynamic stratification 
of smaller scale, more realistic numerical models where the frequency-dependent 
radiative transfer equation in LTE is solved in detail. This technique retains the 
computational efficiency of earlier parameterized methods, but does not require 
continual calibration against more realistic simulations. We find that with the 
new method we are able to initiate and sustain a stable convection pattern with 
a distribution of cell sizes and turnover times characteristic of solar granulation. 

The method presents a middle ground between realistic radiative MHD mod- 
els that solve the transfer equation in detail, and idealized models that simply 
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impose a thermodynamic stratification, or ignore the physics of radiative trans- 
port entirely. The motivation for developing this technique is to make feasible 
physics-based large-scale or global parameter space studies of the interaction 
of active region-scale magnetic fields with the small scale fields associated with 
granular convection in a domain that includes both a convection zone and corona. 

Whether the approximate treatment captures enough of the essential physics 
of the system still remains to be seen. The technique is certainly limited by the 
fact that sideways transport is ignored, and important physics of the chromo- 
sphere has yet to be included in the current models. Even so, we are able to 
generate solar-like convective turbulence in a physically self-consistent way, and 
follow the magnetic evolution of structures that thread the interface between the 
convective interior and corona. 

In particular, we presented two simulations that confirm the existence of a 
"separatrix" in the flow of magnetic energy from the interior to the atmosphere 
(see ISteiner et al.i I2008P , and demonstrate that it is the mechanical energy of 
surface convection driven by strong radiative cooling that is the source of energy 
for this divergent flux of electromagnetic energy. In a quiescent region, in the 
absence of strong concentrations of magnetic flux, our models suggest that it is 
the low photosphere that provides the source of electromagnetic energy to the 
chromosphere and corona, not the sub-surface layers. 
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